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Dynamical processes on networks are currently being considered in different domains of cross- 
disciplinary interest. Reaction-diffusion systems hosted on directed graphs are in particular relevant 
for their widespread applications, from neuroscience, to computer networks and traffic systems. 

Due to the peculiar spectrum of the discrete Laplacian operator, homogeneous fixed points can 
turn unstable, on a directed support, because of the topology of the network, a phenomenon which 
cannot be induced on undirected graphs. A linear analysis can be performed to single out the 
conditions that underly the instability. The complete characterization of the patterns, which are 
eventually attained beyond the linear regime of exponential growth, calls instead for a full non 
linear treatment. By performing a multiple time scale perturbative calculation, we here derive an 
effective equation for the non linear evolution of the amplitude of the most unstable mode, close to 
the threshold of criticality. This is a Stuart-Landau equation whose complex coefficients appear to 
depend on the topological features of the embedding directed graph. The theory proves adequate 
versus simulations, as confirmed by operating with a paradigmatic reaction-diffusion model. 

PACS numbers: 89.75.He 89.75.Kd 89.75.Fb 


Networks are undoubtedly gaining considerable impor¬ 
tance in the modeling of natural and artificial phenom¬ 
ena fli- They define in fact the natural playground 
for a large plethora of problems, that assume a hetero¬ 
geneous support for the connections among constituents. 
In the brain, for instance, neuronal networks provide the 
skeleton for the efficient transport of the electric signal 
(^ . The crowded world of cells in general is shaped 
by veritable routes, the microtubules, that result in an 
intricate cobweb of interlinked paths [l2| . The flow of in¬ 
formation on Internet, and its multifaceted applications, 
heavily rely on the topology of the underlying, global 
and local, network of contacts. Human mobility pat¬ 
terns, with their consequences for transportation design 
and epidemic control, configure, at a plausible level of 
abstraction, as effective graphs, linking different spatial 
locations. 


Reactions occur on each node between species that 
populate the examined system. Individual actors 
(molecules, humans, cars or even bits of information) can 
relocate to distant sites, when exploring the network on 
which they are physically confined. This latter process 
is ruled by diffusion on the heterogeneous, network-like 
support, different avenues of transport being available to 
the microscopic entities, as dictated by the adjacency ma¬ 
trix associated to the hosting graph. The non trivial in¬ 
terplay between reactions and diffusion can instigate the 
emergence of spatially extended motifs 17], [iSj , which re¬ 
flect the inherent ability of the system to spontaneously 
self-organize and consequently perform dedicated tasks. 
In general, when space reduces to a regular lattice or a 
symmetric graph, the dynamics is uniquely responsible 


for the onset of the instability which eventually material¬ 
izes in the observed macroscopic and collective patterns. 
These are, for instance, the celebrated Turing patterns 
that, in recent years, have received much attention also 
in light of their applicability on networks iaQ. 


In applications, however, networks are not always sym¬ 
metric, or, undirected, as customarily termed. Often a 
connection between adjacent nodes imposes a specific di¬ 
rection to the journey, thus resulting in a so called di¬ 
rected graph. 


The map of neural connection is manifestly asymmet¬ 
ric, because of the neurons’ physiology Q- In connectome 
models in fact the coarse-grained maps of the brain re¬ 
veal an asymmetric arrangements of connections at differ¬ 
ent spatial scales [iil,[i3. Cytoskeletal molecular motors 
move unidirectionally along an oriented polymer tracks. 
The cyberworld is also characterized by an asymmetric 
routing of the links As traffic is concerned, several 
routes can be crossed in one direction only, thus breaking 
the symmetry between pairs of nodes. When reaction- 
diffusion systems are considered on directed networks, 
topology does matter. Surprisingly, patterns can rise on 
a directed support, even if they are formally impeded on 
a regular, continuum or discrete, spatial medium. The 
mathematics of this process has been recently investi¬ 
gated in Q, where the conditions for the instability are 
obtained in the framework of a standard linear analysis 
calculation. The patterns which manifest as a byproduct 
of the aforementioned instability, reflect however the non- 
linearities which are accommodated for into the model 
and that are, by definition, omitted in the linear analysis 
theory. In other words, the conditions for the emergence 
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of topology driven patterns for a reaction diffusion system 
on a directed graph can be singled out, but the charac¬ 
terization of the subsequent non linear stage of evolution 
solely relies on numerical methods. 

In this paper we aim at filling this gap, by analyti¬ 
cally deriving an effective equation for the evolution of 
the amplitude of the unstable mode, near the threshold 
of criticality. The spatial characteristics of the generated 
patterns owe to the spectrum of the Laplacian operator 
that governs the diffusion process. The analysis builds 
on a multiple time scale treatment originally devised in 
[MEl, and recently reconsidered for the rather spe¬ 
cific case of a reaction-diffusion system placed on top 
of a symmetric network and subject to weak couplings 
0. At variance, we here focus on the case of a directed 
graph and assume arbitrary large diffusion coefficients. 
This latter condition results in a complexification of the 
analytical procedure: the linear calculation is carried out 
in a iV dimensional space, N being the number of nodes 
in the graphs. The extension to non-linear orders proves 
consequently more demanding. A Stuart-Landau (SL) 
equation is eventually derived for the amplitude of the 
unstable mode. Unprecedently, the coefficients of the SL 
equation reflect the topology of the network, the factual 
drive to the instability. Simulations performed for the 
Brusselator model, a reaction-diffusion system of peda¬ 
gogical relevance, confirm the predictive adequacy of the 
analytical solution, obtained in the framework of the ef¬ 
fective SL scenario. 


amined network. More explicitly, = Aij —Sijki where 
ki = Aij represents the degree of node i. and 
Dy are the diffusion coefficients. To make contact with 
the analysis carried out in Q, we shall deal with per¬ 
fectly balanced networks, namely graphs characterized 
by an identical number of ingoing and outgoing links. We 
will then assume that the equations o admit a homo¬ 
geneous stable equilibrium identified as {x*,y*). To save 
notations, it is convenient to define a vector which con¬ 
tains the concentrations Xi and yi at any node location 
i = l,..,iV, namely x = (cci, ..., xat, yi, ..., Con¬ 

sequently, X* will refer to the aforesaid steady state. We 
are here interested in the conditions that yield a desta¬ 
bilization of the homogeneous stationary stable solution 
X*, as follows the injection of a tiny perturbation which 
activates non trivial interferences between diffusion and 
reaction terms. As anticipated above, the directed spa¬ 
tial support matters: it can actively seed an instability, 
which is instead prevented to occur when the problem 
is formulated on a symmetric spatial backing. In the 
following, we shall briefly recall the main steps of the 
linear analysis theory: these are in fact propedeutic to 
the forthcoming developments, which aim at the full non 
linear picture. 


Linear stability analysis 


RESULTS 


Consider a directed network composed of N nodes. 
The topological structure of the network is encoded in 
the asymmetric adjacency matrix, here denoted by A. 
The element Aij is equal to 1, if nodes i and j are con¬ 
nected, or 0 otherwise. Each node i is populated by two 
species, whose concentrations are respectively labeled Xi 
and yi- The species may react or diffuse throughout the 
network, as specified by the following general set of equa¬ 
tions 


d 


dt 



N 

f{xi,yi,fi) -b 

N 

g{xi,y^,ti) + Dy'^A.jyj 


( 1 ) 


where /(•,•,/.*) and are nonlinear functions of 

the concentration, which descend from the specific reac¬ 
tions being at play, /x is a vector of arbitrary dimension, 
where we imagine stored the scalar parameters, as e.g. 
the rates associated to the reactions chain, which appear 
to modulate the process of mutual and self-interaction. 
A stands for the Laplacian matrix associated to the ex- 


Introduce a small inhomogeneous perturbation, Sxi 
and Syi, to the uniform steady state. In formulae, 
{xi, yi) = {x*,y*) + {Sxi, Syi) for i = 1,..., N. Substitute 
the latter ansatz into equations O: Taylor expanding the 
obtained system and packing 5xi and 6yi into the column 
vector u = (i5xi,..., Jxat, (5yi, ..., one immedi¬ 

ately finds the following equation for the time evolution 
of u: 



(L + D)u -f AIuu + Afuuu 


( 2 ) 


where L and D are two 2N x 2N block matrices 


/ /.(x*)Iy 

1=1 

* 

X 

D = 


On ^ 





\ ya;(x*)Ijv 



\ Oy 

DyAj 


with ly and Oy denoting, respectively, the identity ma¬ 
trix and the null matrix of size N. AIuu and Afuuu are 
symbolic notations, mutuated from Q. These are vectors 
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whose ith components respectively read 


(TWuu), = ^ < 




(AAuuu)^ = - < 


j,k£{i,i-\-N} 

E 

j,k£{i,i—N} 

/(x*) 

E — 7 ^—UjUkUi iov i ^ N 

OXiOXkOXi 
j,k,le{i,i+N} ^ 


E 


j,k,lG{i,i—N} 


5^g(x*) 

dxjdxkdxi 


■UjUkUi for i > N 


The study of the stability of {x*,y*) relies on the linear 
part of equation 


(L + D) u = Au 


(3) 


To solve the above linear system, one needs to intro¬ 
duce the eigenvalues and eigenvectors 0^“^ of the 
Laplacian operator iQ. These are solutions of the 
eigenvalue problem for a = 1,..., A^. 

Importantly, when the hosting network is directed, the 
eigenvalues of the Laplacian are complex. This latter 
property is ultimately responsible for the peculiar be¬ 
havior of reaction-diffusion systems placed on asymmet¬ 
ric graphs, as compared to their undirected homologues. 
The inhomogeneous perturbations Sxi and Syi can be ex¬ 
panded as: 

N N 

Sx, = Sy, = ^ V? (4) 

Oi—l Q;=l 


stable patterns develop when the instability takes place 
on ordinary continuum space or on a symmetric graph. 
These are the celebrated Turing patterns, that typify on 
networks as a material segregation in activator rich and 
activator poor groups. For reaction-diffusion systems on 
directed supports, waves are instead obtained as the late 
time echo of the instability. 

Starting from these premises, we here wish to address 
the full non linear dynamics that stems for a topology 
driven instability, and eventually obtain a close form so¬ 
lution for the emerging traveling waves. To reach this 
goal we shall initialize the system right at the thresh¬ 
old of the instability (/i = /Xq), when the real part of 
the dispersion relation is about to cross the horizontal 
axis, and then perturb the reaction parameter /Xq so as 
to make the homogeneous fixed point slightly unstable. 
A multiple time scale perturbative analysis, which ac¬ 
commodates for key topological ingredients, will open up 
the avenue to a detailed characterization of the complete 
non linear picture. 


When /X = /Xq, the maximum value of is therefore 
identically equal to zero, for a critical index a = ac, 
to which corresponds a selected Laplacian eigenvalue 
A(“A = A^g°^ -I- xAj^^ Since ^ 0, it follows Q 

that A^"”) yf 0. Indeed, A(““i = ±xwo, where wq = 

With A(A^“;)) = {2D,DyK^^:^ + UDy + 


gyDx)/ /a; + gy + {Dy -|- ^ 


as determined from 


a straightforward calculation. From equation ([S]), one 
can readily obtain = -(fx + DxA^^°^)/fy + i{uJo - 
DxA^^)/fy. The solution of the linear problem: 


where Cq depend on initial conditions, and will be 
self-consistently specified later on. By inserting ([4]) into 
©, yields N copy of the following system 


ffx + DxA^<^) 

\ 9x 


fy 

gy + DyA(‘^^ - A(“), 



which admits a non trivial solution provided 


f fx + - A(“) 


det 


= 0 


9y 


+ - A(“) 


( 6 ) 

Equation ([6]) returns a second order polynomial for 
A*^“^ as a function of A*^“\ known as the dispersion re¬ 
lation. The stability of {x*,y*) depends on the sign of 
the real part of A(“\ here termed A^j: if A^^ is negative 
Va, the {x*,y*) is stable, while it turns unstable if A^^ 
crosses punctually the x-axis. In this case, the imposed 
perturbation grows exponentially, in the linear regime of 
the evolution, and the system displays self-organized pat¬ 
terns at the non-linear stage of the evolution. Stationary 


— l 2 Ar — L — D ) u — 0 


(7) 


is hence given by 


u = Uoe““‘ -h c.c. (8) 

where c.c. stands for the complex conjugate. Here Uq = 
is the right eigenvector of L -|- D 
corresponding to the eigenvalue iojQ. As we shall see, 
Uq encodes the spatial characteristics of the predicted 
pattern. 


Multiscale analysis: a topology dependent Stuart 
Landau equation 

Let us start from the neutral condition highlighted 
above, when the parameters are set to the marginal value 
/Xq that yields A^g°^ = 0. Imagine to impose an appro¬ 
priate perturbation in the form /x = /Xq -|- e^/Xj^, where e 
plays the role of a small parameter, and Hi is order one. 
This modulation endows a tiny instability to develop: the 
dispersion relation acquires therefore a positive real part. 
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which consistently scales as e^. This latter observation 
sets the characteristic time scale for the examined insta¬ 
bility, and opens up the perspective for a formal math¬ 
ematical investigation. Following the prescription of the 
multiple time scale technique, we introduce r = the 
slow time variable, which we treat as independent from 
time t. In the solution of the perturbation problem, the 
additional freedom introduced by the new independent 
time variable will be exploited to remove undesired sec¬ 
ular terms. As we shall see, the latter set constraints 
on the approximate solution, which are called solvability 
conditions. 

The total derivative with respect to the original time 
t rewrites: 


d _ d 2 '9 

dt ^ dt ^ dr 


(9) 


Moreover, one may assume the following expansions to 
hold 


L — Lq -|- e^Li -|-... 

M= Mo + €^Mi + ... (10) 

JV= iVo + + ■■■ 


the unperturbed parameters /Xq, and the associated cor¬ 
rection factors fii, being implicitly contained in the def¬ 
inition of the above operators. We further assume that 
u, the solution of the non linear equation can be ex¬ 
pressed as a perturbative series, function of both t and 
r: 


u(t) = eui(t,r)-H e^U 2 (t, r)-H ... ( 11 ) 

To proceed in the analysis, one inserts equations ®, 
m and CD into @ to get: 

/ d d \ 

( ^^2Ar + — Lq — D — e^Li “ • ■ • j (eui-|-e^U2-l-. 

= e^AloUiUi -I- e^(2AloUiU2 -I-A/qUiUiUi) -|- 0{e*) 

Equating terms of the same order in e returns the follow¬ 
ing family of equations 

— Lq — Ml, = IBi, (12) 

with 1 / = 1,2,3.... Following the Fredholm theorem (see 
Appendix), the linear system (1121) admits a non trivial 
solution if the solvability condition is satisfied, namely 
if ((U5)t,BW) = 0, where the angular brackets denotes 
the scalar product. 

We shall hereafter focus on the first three equations 
of the above hierarchy. The corresponding right-hand 
sides (see also Appendix) respectively read Bi = 0, 
B 2 = MqUiUi and B 3 = (-^l 2 Ar-l-Li)ui + 2 A 4 oUiU 2 -l- 
A/qUiUiUi. The solvability condition is naturally met for 
12 = 1 , 2 , while it needs to be explicitly imposed for v = S. 


Consider first the leading order contribution, v = 1 
and solve the corresponding differential equation for Ui. 
As expected, this is equivalent to equation that we 
derived under the linear approximation. Hence, Ui fol¬ 
lows from (|S]) modified with the inclusion of an arbitrary, 
complex and so far undermined, amplitude factor W(t), 
function of the slow time scale r. In formulae: 

ui(t, r) = lF(r)u = lF(r)Uoe*‘^“‘ + c.c. (13) 

As we will see, the factor W{t) sets the typical ampli¬ 
tude of the emerging patterns: it should be constrained 
to match the required solvability condition and so self- 
consistently determined via the multiple scale calcula¬ 
tion. As already emphasized, equation (TT^ constitutes 
a natural generalization of the linear solution ([5]) , which 
indirectly accommodates for the non-linearities through 
the slow varying amplitude factor W. This will in turn 
enable us to track the time evolution of the patterns, be¬ 
yond the initial stage of the exponential growth. The 
remaining part of the calculation is devoted to deriving a 
consistent equation for the time evolution of the complex 
amplitude W. As we shall see, this amounts to imposing 
the solvability condition at = 3. 

To solve the next-to-leading order (z/ = ^ equation in 
(HD, we put forward the following ansatz [Sj for U 2 : 

U 2 = V+e 2 “'>‘ -h -h noui -H Vq 

The constant vq cannot be determined at this stage, and 
will not affect the forthcoming developments. Inserting 
in (fT^ and grouping together the terms that do not de¬ 
pend on t, one finds Vq = -2\W\^ (Lq -b D)"^ TWoUqUo 
where the bar stands for the conjugate. Similarly, equat¬ 
ing the terms proportional to (resp. yields 

V+ = V_ = (2iuJolN -Lo- D)-^ MoUoUq. 

.) At the next order in the hierarchy, 1 / = S, the linear 
equation for U 3 builds on the above characterization for 
both Ui and U 2 . In particular, the unknown complex am¬ 
plitude W enters the definition of the right-hand side B 3 . 
By imposing the solvability condition ((Uq)'I', Bg^^) = 0, 
and carrying out a straightforward manipulation, one 
eventually obtains the following SL equation for W (r) 

JLw{t) = (14) 

dr 

where = g^^l + igfl^ = (Uo)'l'LiUo and 

3(1) ^ = -(U 5 )t [2M0V+U0 

-|-2AloVoUo + 3 A/ 0 U 0 U 0 U 0 ] are complex num¬ 
bers. Notice that g^^^ and g^^'> depend both on the 
reaction terms of the original system o, through e.g. 
Li, A4o, A/q, and on the topological characteristics of 
the embedding network, via Uq. To derive equation 
da use has been made of the normalization condition 

(us)tus = 1 . 
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The solution of equation (fl3 can be cast in the form: 


W{t) = 


(0) 

9Re 


\ \9 


( 1 ), 
Re I 


exp 




. 

where is a phase term which relates to the assigned 
initial conditions. 

Summing up, and recalling equation (fT^ . the wave¬ 
like pattern x(t, e), instigated by the directed network, 
close to the threshold of instability, will be approximately 
described by: 


x(t, e) = x*-|- 


M ig2.’i 


Uq exp 


iuiot ' 



-I- c.c. 


where we have arbitrarily set tp = 0. As anticipated, the 
structure of the graph which ultimately drives the insta¬ 
bility enters parametrically the above solution (1161) . In 
the following, we shall indicate with the amplitude 
of the oscillating patterns, for respectively species x and 
y, around the average solution. 


Alternative perturbation scheme: acting on the 
diffusion coefficients 

In the previous section we have seen how to charac¬ 
terize the emerging patterns when they originate from 
a perturbation of the reaction coefficients fi. Similarly, 
one could imagine to induce the instability by perturbing 
the diffusion constants from and Dy. More specifi¬ 
cally, we initialize the unperturbed system so as to match 
the marginal condition = 0 and then perform the 

change ^ and Dy ^ Dy + e^Dyi, where 

e is a small parameter, and D^i and Dyi are order one 
scalar quantities. Proceeding in analogy with the above 
yields a SL differential equation for the evolution of the 
complex amplitude factor W(t), where gi is unchanged 
and go = (Uq)'*'DiUo where 


( D,iA 

Oat ^ 

\ Oat 

DylAj 


Numerical validation of the theory 

We here aim at testing the predictions of the theory, by 
drawing a comparison with the outcome of direct simula¬ 
tions performed for a reaction-diffusion model of paradig¬ 
matic interest. This is the celebrated Brusselator model. 


a non linear reaction scheme which describes the auto- 
catalytic coupling of two mutually interacting chemical 
species. Details of this model can be found in the Ap¬ 
pendix. 

As a first example, we consider the Brusselator model 
defined on a balanced network generated with a slightly 
modified version of the Newman-Watts (NM) algorithm 
[l6| (see Appendix). In the left panel of Figure [T] we 
display with symbols the real part of the dispersion rela¬ 
tion ^ function of the real part of the Laplacian 

eigenvalue A^^ (changed in sign). The parameters of 
the model have been set so as to have the largest value of 
equal to zero, in correspondence of a specific 

) The solid line represents instead the dispersion relation 
obtained, with the same choice of the parameters, for the 
limiting case of a symmetric continuous support. If the 
(lS)^stem is placed on top of a symmetric graph, the con¬ 
tinuous curve turns into a discrete collection of points, 
following exactly the same profile and reflecting the finite 
set of (real) eigenvalues, associated to the Laplacian oper¬ 
ator. When the embedding network is instead asymmet¬ 
ric, the complex component of the Laplacian spectrum 
lifts the dispersion relation, as depicted in leftmost panel 
of Figure [U so eventually inducing a topology driven in¬ 
stability, in a otherwise dynamically stable system. In the 
other two panels of Figure [1] the patterns obtained via a 
numerical integration of the reaction-diffusion system (HD 
and the analytical solution (ITCl) are respectively reported, 
displaying a satisfying degree of correspondence. 

As an additional check for the developed theory, we 
consider a family of directed regular lattices, with vary¬ 
ing level of imposed asymmetry. More specifically, we 
preliminary assumed a closed one dimensional ring com¬ 
posed of N nodes: each node has K links to its first K 
nearest neighbors encountered when circulating the ring 
clockwise. The adjacency matrix which describes such 
a lattice is then shifted, via n successive applications of 
a one-step shift operator, so to result in a set of dis¬ 
tinct lattices, which tend to progressively approach the 
symmetric limiting case. For such n directed networks, 
we computed the amplitude of the predicted, topology 
driven patterns, as follows equation dH, and compared 
it to the outcome of numerical simulations based on the 
original reaction-diffusion model. Results of the analy¬ 
sis are reported in Figure [U where the amplitude of the 
pattern is plotted as a function of the degree of shift n. 
Here, the instability is produced upon perturbation of 
the reaction parameter b. An overall excellent agreement 
is observed, between theory and simulations. The pre¬ 
dictive adequacy of the theory can be also appreciated in 
Figure[3]where the time dependent patterns are displayed 
for n = 5, a representative case study. The same conclu¬ 
sion holds when the perturbation acts on the diffusion 
coefficient D^ and Dy (data not shown). 
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DISCUSSION 

Self-organized patterns can spontaneously develop in a 
multi-species reaction-diffusion system, as follow a sym¬ 
metry breaking instability of an homogeneous equilib¬ 
rium. Inhomogeneous perturbation can in fact amplify 
due to the constructive interference between reaction and 
diffusion terms, and eventually yield coherent, spatially 
extended motifs in the non-linear regime of the evolution. 
Reaction-diffusion systems placed on symmetric graphs 
have been also analyzed in the literature. The conditions 
for the deterministic instability are derived via a linear 
stability analysis, which requires expanding the pertur¬ 
bation on a complete basis formed by the eigenvectors 
of the discrete Laplacian. For system hosted on undi¬ 
rected networks, the instability is essentially driven by 
nonlinearities, which stem from both reactions and diffu¬ 
sion. The topology of the embedding network-like sup¬ 
port defines the relevant directions for the spreading of 
the perturbation, but cannot influence the onset of the 
instability. A radically different scenario is encountered 
when a directed graph is instead assumed to provide the 
spatial backing for the scrutinized model. In this case, 
the topology of the space is equally important and sig¬ 
nificantly impact the conditions that drive the dynamical 
instability. 

Building on these recent advances, the aim of this pa¬ 
per is to go beyond the standard linear stability analysis 
for the outbreak of the instability and provide a com¬ 
plete characterization of the patterns emerging on a di¬ 
rected discrete support, in the fully developed non linear 
regime. To this end we have applied a multiple time scale 
analysis, generalizing to the present context the original 
derivation of [I3- This results in a cumbersome calcula¬ 
tion owing to the particular nature of the diffusive cou¬ 
pling imposed. The amplitude of the most unstable mode 
is shown to obey a Stuart-Landau (SL) equation whose 
coefficients unprecedently reflect the topology of the net¬ 
work, the genuine drive to the instability. Simulations 
performed for the Brussellator model, confirm the valid¬ 
ity of the theory, which proves effective in quantitatively 
grasping the characteristics of self-emerging dynamical 
patterns, close to the threshold of instability. This is 
a significant achievement which could translate in novel 
strategies to control the dynamics of the system, via ap¬ 
propriate fixing of topological features, including the su¬ 
pervised addition/removal of specific nodes/links in the 
network. 

APPENDIX 
The solvability condition 

Let A be a linear operator, and u(t) and b(t) two 
complex vectors of the same length. According to the 


Fredholm theorem, a linear system Au(t) = b(t) is solv¬ 
able if {v(t),b(t)) = 0 for all vectors v(t) solution of 
A*v(t) = 0, where A* is the adjoint operator satisfy¬ 
ing (A*y,x) = (y. Ax) 'ix,y. The angular brackets de¬ 
note the scalar product that we here define as (v(t), b(t)) 
= vl(t)b(t)dt, the symbol f standing for the con¬ 

jugate transpose. With reference to equation ([12]), the 
first requirement of the Fredholm theorem consists in 
finding v(t) such that {d/dtl 2 N — Lq — D)*v(t) = 0. 
Recalling that Lq and D are real matrices, by par¬ 
tial integration we find that {d/dt\ 2 N — Lq — D)* = 
-{d/dtl2N + Lo + D)^. As a consequence, the sys¬ 
tem to be solved is —{d/dtl 2 N -I- Lq -I- D)’^v(t) = 0. 
In analogy with equation (0, we search v(t) in the 
form v(t) = UQe“°‘ for some vector Ug. Substituting 
this ansatz into the previous equation, we find (Lg -I- 
D)^Ug = —zwoUg. In analogy with Uq, Uq is related to 
the eigenvalue problem = (A^^ — 

through Uq = with = 

-ifx + D^A^^^^)/gy-i{u}o-D^A\°^^)/gy. Having dehned 
US, one can explicitly write the solvability condition 
(USe®‘^“*, B^(t, t)) = 0. Since By(t, r) turns out to be pe¬ 
riodic functions of period 27r/a;Q, it is appropriate to ex¬ 
press them in the form B^(t,r) = Et -00 
If we multiply this series by (USe“”*)^ we again obtain 
periodic functions that, when integrated over the period 
27r/a;o give zero. The only exception holds for I = 1 which 

gives (USe“'>*,B[,^V)e*“”‘) = (U*o)^Bi^\T)dt. 

The integrand does not depend on time t and therefore 
the integral is zero only if the integrand itself is identi¬ 
cally equal to zero. For this reason the solvability condi¬ 
tion reduces to (US)^b[,^^(t) = 0 Viz. 

The Brusselator model 

In the Brusselator model, the two reaction terms are 
specified by f{xi,yi,g) = a — {b + d)xi -\- cx^yi and 
g{xi,yi,fi) = bxi — cx^yi, where g = (a, 6 , c, d) defines 
a set of positive real parameters. The unique homoge¬ 
neous equilibrium point is {x*,y*) = {a/d,bd/c/a). 

Network generation strategy 

We start from a substrate AT-regular ring made of N 
nodes. The NW algorithm 0 is designed to add, on av¬ 
erage, NKp long-range directed links, in addition to the 
links that originate from the underlying regular lattice. 
Here p lies in the interval [0,1] and represents a proba¬ 
bility to be chosen by the user. The NW algorithm here 
employed is modified so as to result in a balanced net¬ 
work (identical number of incoming and outgoing links, 
per node). To this end, the inclusion of a long-range link 
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starting from node i is accompanied by the insertion of 
a fixed number (3 is our arbitrary choice) of additional 
links to form a loop that closes on i. 
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FIG. 1: Panel (a): The (blue, online) triangles represent the real part of the dispersion relation as a function of (minus) the 
real part of the eigenvalues of the hosting directed network. The black line originates from the continuous theory. Panel (b): 
Pattern emerging from species y as obtained by direct integration of system ©. The concentration on each node is plotted as 
a function of time. Panel (c): Pattern relative to species y as determined from the analytical solution of (I16II . The network is 
made of = 100 nodes and has been generated following the NW recipe with p = 0.27. Parameters are a = 2.1, = 4.002, 

c = 1, d = 1, Dx = 1 and Dy = 3. The perturbation is here acting on & as 6 = &c(l + e^), with = 0.1. 



FIG. 2: Amplitude of the self-emerging oscillations as a func¬ 
tion of the shift parameter n. Here K — 27 and n € (0,12). 
Orange dots refer to the amplitudes obtained from equation 
(nsj, while green symbols follow from numerical integration 
of the Brusselator model. The main panel refers to the ampli¬ 
tude of the patterns relative to species y, while, in the inset, 
the amplitudes are calculated for species x. Parameters are 
a = 4, c = 1, d = 1, Dj: = 1, Dy = 3. The instabilities come 
from the perturbation of parameter b as b = bc{l + e^) for 
= 0.1. For each n, the critical value of be is calculated so 
as to satisfy the condition max{\^^l) = = 0 . 
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FIG. 3: Waves on a directed lattice obtained by imposing n = 5 shifts. For the parameters’ description refer to the caption of 
Figure (2] Panels (a) and (b) show the time dependent patterns relative to species x obtained, respectively, from the original 
Brusselator model and the SL equation. In panel (c) we display x versus the index of the node, while the inset reports a; as a 
function of time t, for a selected node. In both cases, green circles refer to data extracted from panel (a) (simulations), orange 
squares to panel (b) (theory). 
















